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This paper deals with the solution of the spherically symmetric time-dependent Hartree-Fock 
approximation applied in the case of nuclear giant monopole resonances in the small and large 
amplitude regimes. The problem is spatially unbounded as the resonance state is in the continuum. 
The practical requirement to perform the calculation in a finite-sized spatial region results in a 
difficulty with the spatial boundary conditions. Here we propose a absorbing boundary condition 
scheme to handle the conflict. The derivation, via a Laplace transform method, and implementation 
is described. The accuracy and efficiency of the scheme is tested and the results presented to support 
the case that they are a effective way of handling the artificial boundary. 
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I. INTRODUCTION 

It occurs in many areas of physics that the time- 
evolution of a spatially unbounded system is required to 
be analysed. Such systems have been studied in many 
fields of physics involving wave propagation, spanning 
areas such as laser physics and gravitational waves, [1- 
4]. Examples occur in nuclear physics and we analyse 
such a case in the present work. 

The particular physical phenomenon being studied 
here is the nuclear giant monopole resonance. It is well 
known that these are above the particle decay threshold 
[5], so that one allowed decay mode involves the expul- 
sion of one or more nucleons from the nucleus. A time- 
dependent simulation of such a decay will involve the 
spatial region in which the nuclear wavefunction is non- 
negligible becoming larger and larger as time goes on. 

One way of analysing this sort of system is via the time- 
dependent Hartree-Fock (TDHF) approximation, that 
reduces the many-body interaction to a simpler mean 
field one. The simplification however still does not allow 
analytic solutions to be gained but allows for numerical 
analysis to be applied and the computational cost to be 
manageable. 

A common numerical implementation is to discretise 
the equations using time and space grids employing finite 
difference methods. Here, a non-trivial problem occurs 
because the boundary of the finite grids impose an arti- 
ficial boundary into the solution. As the outgoing wave 
condition for the Hartree-Fock equation is evaluated at 
infinity it cannot straightforwardly be applied directly. 
Enforcing the wrong boundary conditions results in the 
solution becoming incorrect for the time after the emitted 
particles have reached the artificial boundary and so it 
can be important that the boundary is handled properly. 

There are various methods available that aim to sim- 
ulate or circumvent the application of the outgoing wave 
condition [6, 7]. The most crude is just to apply a reflect- 
ing boundary sufficiently far away so that the matter 
being emitted does not reach it within the time of the 
calculation. This works and reflecting boundaries can be 
easily implemented but the major drawback is that one 
needs an increasing number of grid points in space as one 



wants to evolve further in time. Eventually, this becomes 
computationally unfeasible. 

Other methods include absorbing potentials and mask- 
ing functions. These allow the artificial boundary to be 
placed closer to the nucleus but generally have to be 
tuned to each particular case and do not in general ap- 
proximate the outgoing wave condition perfectly. 

Here, we present a method of implementing exact 
boundary conditions [1, 2]. These rely on choosing the 
artificial boundary such that the potential outside of it 
has a simple form, so that the propogation of waves in the 
exterior region does not have to be dealt with explicitly. 

In solving the TDHF equations, a simplified Skyrme 
interaction is used in the implementation which repro- 
duces the magic numbers needed for ^He, ^O and ^Ca 
to be seen without the complexity of the full interac- 
tion [8], as a reasonable proof-of-concept. Spherical sym- 
metry is also assumed inside and outside of the artificial 
boundary. The calculations involves one, in the case of 
JjHe, or more, in the cases of ^O and ^[jCa, different 
forms of differential equation, each of which requires its 
own absorbing boundary condition to be applied. Here 
some continuous absorbing boundary conditions are used. 
Other types of absorbing boundary are fully-discrete [9] 
and semi-discrete [10] but are not described here. A re- 
view of the various absorbing boundary conditions can 
be found in [11]. 

The structure of this paper is as follows: Section II 
gives a brief summary of nuclear giant monopole res- 
onances; sections III and IV describe the Hartree-Fock 
approximation, the first the theory and the second its 
discretization and implementation; section V and VI de- 
scribe the exterior problem and the absorbing boundary 
conditions; sections VII, VIII and IX show the testing 
and results of our implementation which includes a short 
analysis of the errors caused by the discretization and 
strength functions for *He, ^O and ^[jCa, and results 
with large-amplitude excitation. We end with some con- 
cluding remarks. 
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II. GIANT MONOPOLE RESONANCES 

Giant resonances are collective modes of excitation of 
finite fermionic systems [12]. The first evidence for their 
existence in atomic nuclei came in 1937, with a the- 
oretical description and systematic experimental study 
coming in the next decade [13]. While the first studies 
excited the electric isovector dipole resonance, in which 
protons and neutrons oscillated out of phase with each 
other due to the dominance of the El component of 
the photon field, other giant resonances were discovered 
later. In particular, the isoscalar giant monopole reson- 
ance (GMR) was definitively reported in 1977 [14]. 

The GMR, as a compression mode, probes the nuc- 
lear equation of state [15], and is therefore useful in con- 
straining nuclear models [16]. As a spherically-symmetric 
excitation, it is the first port of call for testing new the- 
oretical methods, as the symmetry renders many types 
of calculation more simple. In particular, methods based 
on Time-Dependent Hartree-Fock have turned to giant 
monopole resonances in spherical doubly-magic nuclei as 
a proving ground [17-23]. 

The present paper is written in that spirit, employ- 
ing the simplified t^-t^ version of the Skyrme force used 
in previous applications [18, 20]. While the focus of 
this work is on the development of the boundary con- 
ditions, and the simplified Skyrme force we use should 
not be expected to give good agreement with experiment, 
it is noted that of the three nuclei considered here, the 
GMR has been unambiguosly observed only in 40 Ca [24], 
though the nature of giant resonances in general in nuclei 
as light as 4 He is a subject of ongoing interest [25]. 

The key observable calculated for the giant resonance 
is the linear response function, describing the response 
of the nucleus to an external perturbation [26]. From 
this, one derives the strength function, related in turn 
to the experimental cross section for the reaction. The 
strength function can be obtained, within TDHF, via the 
Fourier Transform of the time-dependent moment of the 
resonance mode desired [27] and we present calculations 
of such strength functions. We note that the strength 
functions are particular sensitive to the success of imple- 
mentation of the absorbing boundary conditions [6] , and 
provide a good measure of success, as well as being the 
physically relevant quantity. 



III. TIME DEPENDENT HARTREE-FOCK 

The time-dependent Hartree-Fock method originates 
with Dirac [28], and became computationally viable for 
nuclear processes in the 1970s [29-31]. Since then it 
has been extensively used for calculating heavy-ion re- 
actions [32] and giant resonances [33], with increasingly 
sophisticated implementations of the effective interaction 
[34, 35] . A full derivation of the Time-dependent Hartree- 
Fock equations in the case of Skyrme forces can be found 
in the original paper by Engel et al. [36]. In the present 



case, with the simplified Skyrme force, and omitting Cou- 
lomb, we note that the Time-Dependent Hartree-Fock 
equations can be written as a series of coupled non-linear 
Schrodinger equations of the form 



ih = hip\[r,t), 

at 



A = l, 



■ A (i) 



where the Hartree-Fock Hamiltonian is given by 

h 2 

h = - — V 2 + ap(r,t) + bp 2 (r,t) 1 (2) 
2m 

with p(r, t) = X^a=i ^a(^ *)^a(^) denoting the particle 
density. The values of a and b used thoughtout this paper 
are taken from [37] where they have the values —817.5 
MeV fm 3 and 3241.5 MeV fm 6 . In practice, the time- 
dependent Hartree-Fock equations are solved by evolving 
in time according to 



(3) 



Specialisation to spherical symmetry, and details of dis- 
critisation methods, are given in the following sections, 
in which the details of the algorithm dealing with the 
boundary conditions are also given. 



IV. INTERIOR DISCRETIZATION 

As well as the coupled non-linear differential equations 
noted in the previous section, initial conditions are re- 
quired, and are calculated from stationary Hartree-Fock. 
We first describe our method for calculating the station- 
ary solution and then go on to the time-dependent case. 
In both we discetized the equations on equally spaced 
grids, for simplicity, though non-uniform grids can be in 
themselves useful in pushing the boundary far into the 
exterior region at an acceptable computational cost [38]. 



A. Stationary Discretisation 

We start with the calculation of the initial condition, 
which itself is a non-linear problem. We solve it by the 
following iterative procedure: 

M 4) (r)Qi l+1) (r) = A£ i+1) QL 4+1) (r) 



1 j9f 
"2 dr- 



i Q (io + l) 
2r 2 



+ V{p«(r)} 



p {0 Kr) = I ^E a g a \Q i " ) (r)\ 



(4) 
(5) 
(6) 

for i e No and where Q(r) = rip represents the reduced 
wave function, V is the potential, and l a the orbital angu- 
lar momentum. We calculate the initial guess, p^°'{r, t), 

using harmonic oscillator wave-functions as the Q^P (r) 
in equation (6). 

Spatial discritisation of the equations is made on a 
uniformly-spaced grid, such that 

Rout 



mAr, m — 1, ... , M, Ar 



M 



(7) 
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where M is the total number of gridpoints and R ou t is the 
distance from the origin to the spherical outer boundary. 
The second derivative operator in (5) is treated with the 
three-point approximation. 

We also require the wave functions at two addi- 
tional points; Qa\r ) = Qa (0) and Qa {tm+i) = 
Qa? ((M + l)Ar). Although our differential equation is 
not evaluated at these points, values of the wave function 
here are needed for the finite differencing. 

Working with the reduced wave function leads to a 
boundary condition of Qa\r ) = 0. However, the large- 
r boundary condition, that the wave function remain 
square-integrable and fall to zero strictly only at infin- 
ity and cannot be applied directly. We make use of that 
property that the wavefunctions for bound states decay 
exponentially as r increases. Hence we can find a radius 
at which the wavef unction is zero, within a given accur- 
acy, and so we choose Qa (fM+i) = for the solution of 
the static Hartree-Fock equations. 

This leaves us with a tridiagonal matrix eigenvalue 
problem at each iteration, which can be solved efficiently 
using the LAPACK subroutines. 

We iterate until both the eigenvalue, A« , and the 
mean square errors for each wave function, 



(Q^ +1) \H {i) \Q { l 



-<q£ +1) i(# w ) 2 |q2 +1) > 



(8) 



have stopped changing, within machine precision, from 
one iteration to the next. 



B. Time-Dependent Discretisation 

After the initial states have been found using the above 
procedure we need to apply the monopole boost operator 
in order to start the nucleus in the breathing mode. This 
can be done using the usual boost operator for an iso- 
scalar monopole mode 



Qa (j*rn) t 



(9) 



where k is the adjustable strength. 

Once this has been done the QaS can be propagated 
in time. The equally spaced time grid 



nAt, n = 1, 



,N 



(10) 



is used and the same space grid, (7), as the stationary 
problem. The Crank-Nicholson method is then used for 
the time discretization of the time-dependent Hartree- 
Fock equation: 

iAt 



(11) 



I + -^-H(r m ,t n _i) Q a (r,t n ) 



iAt 



H(r m ,t n _i) Q a {r,t n -i) 



We choose the Crank-Nicholson method because it has 
properties that are useful for this type of calculation: it is 
unconditionally stable; and it maintains norm. However 
being an implicit method it also yields the Hamiltonian 
evaluated at a half time-step and so through the potential 
term the density evaluated at the half timestep. This 
means our resulting equations are not a system of linear 
equations. To get around this problem we use an explicit 
method, which is calculated after each propagation in 
time to yield the wavefunctions needed to calculate the 
half-time-step density. We use a method based on the 
evolution operator: 



iAt - 

Q(r m ,t n+ i) = exp ( -—H(r m ,t n ) ) Q(r,t n ) (12) 



2 -f 

j=0 v 

+£>(Ar 2 , At Jma:c ) 



H(r m ,t n )) Q(r m ,i r J13) 
(14) 



requiring knowledge of the Hamiltonian only at the cur- 
rent time-step. 

Once equation (11) has been discretized in space us- 
ing central differences and the grid (7) it is a tridiagonal 
matrix equation, again solved with LAPACK routines to 
get from one time to the next. 

However, the last row in the matrix contains an un- 
known Q(rm + i,t n ) for n > 0. This has to be specified 
with the boundary condition which we know at infinity, 
but we require a boundary condition at r = (M + l)Ar. 
We could use the same reasoning as the stationary case, 
that we can find a point at which the wavefunction will 
be zero and apply the boundary there. We also know 
however that this system has a probability of particle 
emission, which manifests itself in the calculations as a 
thin non-zero tail travelling away from the central mass 
near the origin. This means as time passes the point 
at which the wavefunction is zero gets increasingly fur- 
ther away. This corresponds to longer calculation times 
which can be prohibitive. Hence we seek an absorbing 
boundary condition to give the value of Q(rM+i,t). 



V. PROBLEM IN THE EXTERIOR 

A. Splitting the Domain 

We start by splitting the domain into two regions: an 
interior in which we choose to contain all the nuclear 
dynamics; and an exterior where we assume only the 
long ranged components are of significance, in this case 
just the centrifugal barrier. Given the partial differential 
equation for a single particle state in coordinate space: 



+0(Ar 2 ,Ai 2 



d 

i-^Qi(r,t) 



Id 2 
2d7 



2 +V(r,t))Q l (r,t), 



(15) 
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with boundary conditions: 



Q,(0,t) = 0, 
lim Qi(r,t) = 0. 



(16) 
(17) 



We can mathematically describe the splitting with the 
potential term: 

V(r,t) = V short {r,t) + V long (r), (18) 

where we define: 

V short (r,t)=0ioTr>R, (19) 

(20) 



V long (r) = for r ■ 0. 



2r 2 

The problem has now been split into where the internal 
potential is present and where it is not. The parameter 
R is commonly called the artificial boundary and has to 
be chosen so equations (19) and (20) are satisfied. We 
also assume that the initial wave function is zero outside 
the artificial boundary: 



h{r,0) = for r > R. 



(21) 



This is not overly restrictive and consistent with our 
choice for the solution of the static Hartree-Fock equa- 
tions. 



B. Deriving the Absorbing Boundary Conditions 

We have now all the assumptions needed to construct 
the absorbing boundary condition. There are various 
ways of doing this and a Green's function approached 
has already been described by Heinen and Kull in [1, 2] 
for this problem. We proceed differently, however, by de- 
scribing a derivation using a Laplace transform method. 

We start by recalling the definitions of the Laplace 
transform[39, Chapter 29] in time, f(s), of a function, 
/(*), as: 



f(s) = / f(t)e- st dt, 
Jo 



(22) 



and the inversion formula, known as the Bromwich 
integral [40]: 



/(*) 



1 

2~7ri 



f(s)e st ds. 



(23) 



Combining equations (15) and (20) for r > R we have: 



\&>_ 
'^dr 2 



1(1 + 1) 

2r 2 



Qi(r,t), (24) 



Multiplying by e~ st and integrating time from to oo 
allows us to use equation (22) to get the ordinary differ- 
ential equation: 



l d 2 Qi(r, s) 
2 dr 2 



The substitution Qi(p,s) = phi(p,s) where p = kr and 
k = y/2is, yields the following equation for hi(p, s): 



,d 2 h 



. t>i , r2r-^ + (r*-l{l + l))h l = Q, (26) 

where the square root is assumed to be on the branch 
resulting in a positive real part. As / € No we can see that 
this equation has spherical Bessel functions as solutions 
[39, Chapter 10] of which there are various satisfactory 
pairs. We choose the particular solutions as the spherical 
Bessel functions of the third kind, also known as spherical 
Hankel functions. Any pair of solutions can be used to 
give the same end result once the boundary condition 
are applied. However this pair simplifies the consequent 
derivations. 

Taking the Hankel function solutions, we can write Qi 
as: 



&(r,*) = A(s)ph ( i 1) ( P ) +B(s)phY > {p) 



(2), 



p—kr 



(27) 



Only the boundary condition (17) is relevant here, to 
be precise its Laplace transform, as r > R and may be 
applied by the use of the following limiting forms for 
z — > oo: 



hf\z) 



i- l - l z~ x e iz , 



i'+ 1 z- 1 e—. 



(28) 
(29) 

Assuming c > in the Bromwich integral (23) allows us 
to say that y > where k = \j2is = x + iy, along the 
integration path. So by the limiting form of Q(r,s) as 
r — > oo: 

Ql(r, s) - A(s)r'-V ix -^ r + B{s)i l+1 e {y - tx)r , (30) 

we must have B(s) = 0. 

Q(r,t) and its r derivative can now be written as: 



dQi(r,s) 
dr 



h(r,s) = A(s)ph < i 1 \p) 



p—kr 



p—kr 



(31) 
(32) 



Division of these two equations and evaluating on the 
artificial boundary yields the Laplace transform of the 
absorbing boundary condition: 

/ 



Qi(R,s) 



V 



i_ph ( l 1 \p) 

k d_ 

dp 



K } (p)) 



6Qi(R, s) 
dr 



(33) 



p—kr J 



Use of the convolution theorem for Laplace transforms 
gives us the absorbing boundary condition: 
ft 



Q t (R,t) = g ,(iZ, T ) gg^ T) dr, 



where we define: 

Gi(R,s) 



1 ph^ip) 



k d_ 
dp 



(34) 



(35) 



p—kr 
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G(R,s) being the Laplace transform of G(R,t), which 
can be simplified by the recurrence relation: 



to: 



phf\p) 



Gi(r,s) = 



k{l^l)h^\p)- P h^M 



(36) 



(37) 



p—kr 



C. Calculation of the kernel G(r,t) 

Our final task, before discretization, is to calculate the 
inverse Laplace transform above. This is done by using 
a series expansion[39, p439] for h^\z): 



=i- l - l z-\ 



1 1 

£(Z+-,fc)(-2zz)- fc , (38) 



where: 



1 p + iQl 

V + 2' k) - v\(l-v)V 



(39) 



After manipulation and simplification we gain the ra- 
tional function in k: 



Gi(R,s) 



(i+l 



(;+|,o)(-2ifl)» 



k 



(i+|,0)(-2iR)" 



-{40) 



This can be expanded in partial fractions: 

l+i 



.7 = 1 J 



i + l J 



/2( 



(41) 
(42) 



where the kj are the roots of the polynomial in the de- 
nominator of (40) and a,j are the pole strengths. In prac- 
tice we calculate the roots and strengths for each I with 
Maple. 

The inversion of (42) is performed just by applying the 
well known result from tables[39, 41]: 



v/5 



1 = — = — aw(iaVt), 



(43) 



rather than contour integration of the Bromwich integral 

2 

(23). Here w(z) = e~ z crfc(— iz) is the Faddeeva func- 
tion, which can be calculated with an implementation of 
reference [42]. G(R,s) can now be written as: 



i+i 

Gt(R,T) = J2 

7=1 



'2mt 



-iajkjwizj) 



(44) 



where za 



— -n, 3 y y. Simplification of the above can be 

made by using the limiting form (28) in equation (35) 
and comparing to (41) in the limit k — > oo: 







lim 

k— ^oo 



lim 

k— >oo 



(kGi{r,s)-kGi{r, s) 



P °- p (i-l-^P) 



i+i 

p=kr 3=1 



(45) 
,(46) 



the differentiation of the limiting form is allowed as the 
functions hi (z) are analytic. The limit can be per- 
formed to give: 



i+i 

£' 

i=i 



(47) 



which allows us to write our final form of the kernel G 
as: 



. i+i 



Gj ( i? ' r ) = yp= °A' W ' 



(48) 



An interesting and reassuring feature of this boundary 
condition is that for I = where equation (15) reduces 
to the free one dimensional Schrdinger equation, we have 
the values a± = — i and k\ — 0. Using these values we 
gain the absorbing boundary condition for the free one 
dimensional Schrodinger equation as found in [43]. 



VI. BOUNDARY DISCRETIZATION 

A. Removing the Singularity 

Equations (34) and (48) will now be discretized on the 
grid for use in the Crank-Nicholson scheme. Inspecting 
equation (48) we see that it has a square root singularity 
at r = and is not ideal for numerical integration. So 
integration by-parts is done on the first term to give: 



[2ir d i 1+1 
Gi(R,t) = — — - ~£a^w(^') 

j=i 



(49) 



Our function is now continuous at r = and although 
its derivatives are not it is better suited to the numer- 
ical integration. Note that Gi(R,t) is now an operator. 
Defining a function u^ 1 '(R,t) allows for a more compact 
expression: 



Gi(R,t 



<2iT d 



u {1 \R,t) 



(50) 



7T <9t 
. i + l 

t t W(i? ! r) = -^£a J -%w(^) (51) 



j'=i 
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B. Time Discretization 



Where: 



We first form a semi-discrete equation on the grid t n = 
nAt with t = t/v and r„ = t n . By using the extended 
midpoint rule: 

/ /(r)dr = At^/(t„ +| )+0(At 2 ) (52) 

to evaluate the integral and the difference formulas: 

f{r,t n _i) = + 0(At ) (53) 



/QvU + f(r, t n -i) t /0 _ x 
2 

9/(r,t n _i) _ /(r,t n )-/(r,t„_i) , , ,, 



for functions evaluated at a half time step gives the fol- 
lowing semi-discrete equation: 




7T 2 V ' 5- 



'2^i At ^ dQi(R,t N -i) 



7T 2 ' ' 5' 



f/r 



AT-l 




2ii^A At \ dQi(R,t N - n ) 



2^„+i At s \ dQj(i?, t W -n-i) 



r — + -y^(i?,t n+ i) 



+C(At 2 



C. Space Discretization 

For the space discretization we choose the artificial 
boundary at R = r M _i between the penultimate and 
final spatial grid-points. The following difference formu- 
las are used: 



f(r M ,t) + f(r M -x,t) 



df{r M _x,t) f(r M ,t) - f(r M -i,t) 



0(Ar 2 )(55) 
-C(At 2 )(56) 



dt At 
at the points between the spatial grid. This yields the 
fully discetized absorbing boundary condition: 

(l - B< M ' 0) ) Q,(r Ms %) + (l + B$ M >°^ Q^M-i,*iv) 
= cf M ' 0) (Qi(ru-i,tN-i) - Qi(ru,tN-i) J 

N-l 

^ B\ M ' n) \ Qi{r M ,tN-n) - Ql(rM-l,tN-n) J 



n=l 
JV-1 



53 C\ M ' n) (QiirM-xM-n-x) - Ql — n— 1 ) 

n=l 

0{Ar 2 ,At 2 ). (57) 



.4: 



-2 /iAt 

A7V^T : 



At 



C\ M > n) = .4^2^+1- ^(r„ 4) i„ + .). 

Within the implementation, equation (57) replaces the 
last row of the matrix described in section (IV B). 



VII. RESULTS AND TESTING: ABSORBING 
BOUNDARY EFFECTIVENESS 

Before calculating the giant resonances, the imple- 
mentation of the absorbing boundary is tested in a sim- 
plified case, without any potential, beyond that com- 
ing from the centrifugal term. We apply the absorbing 
boundaries to a partial differential equation of the form 
(15). This is to show the validity of the implementation 
and to demonstrate its performance. The solution to the 
following partial differential equation is found: 



.dQi 



i d 2 Q l (r,t) 

2 dr 2 



1(1 + 1) 
2r 2 



Ql(r,t), (58) 



Qi(r,0) = Are-^l 2 , (59) 
Qi(0,t) = 0, linv-»oo Qi(r,t) = 0, (60) 

for I = 0, 1, 2. Although calculations can be done for any 
angular momentum these are the only values required for 
the Hartree-Fock calculations shown later. A is chosen to 
normalise Qi(r, 0) and is calculated with Simpson's rule. 

Physically the equation corresponds to the evolution 
of a free particle which initially is a shell surrounding 
the origin. Although this sort of system provides no par- 
ticular physical insights, it does allow us to make quick 
and simple calculations which are suitable for testing the 
validity of the method. 

We use the same time and space discretization as de- 
scribed in section (IV) to discretise equation (58). The 
intermediate step (12) is not needed here, as the equation 
is linear. 

Our results will show comparisons between a calcula- 
tion done with absorbing boundaries at r = 10 and one 
with reflecting boundaries at a radius chosen so reflection 
does not occur, which will be specified for each test. 

For this simplified case, we take K = m = 1. 



A. Densities 

To show how the solutions to equation (58) evolve 
through time the probability densities are presented. 
These are gained from calculating the wavefunction 
through time with a reflecting boundary at r = 100. In 
the time interval chosen, [0, 15], reflection does not occur. 



7 



Figure 1 shows us the densities through time for each an- 
gular momentum. Only the interval [0, 10] is plotted as 
this is where we place the test absorbing boundary. The 
results are calculated with grid spacings Ax — At = 0.1. 



In each case we see the bulk of the density begins 
centred at r = 5. As it the system evolves, the wave- 
packet spreads out, and interferes with itself as it reaches 
the origin. 
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Figure 1. These figures show wavefunctions, of angular mo- 
mentum I — changing in time with a percentage leaving 
the interval of interest. The calculations are done with a 
reflecting boundary at r = 100 and have grid spacings of 
Ax = At = 0.1. From top to bottom the graphs show the 
evolution of the wavefunctions at times 0,5,10 and 15. 



B. Radial Comparison of Wavefunction 

We now go on to see how the absorbing boundary per- 
forms. We plot: 



\Q\ ReS \r,t) 



Q\ ABC \r,t)\ 



(61) 



at t = 15, where Q< ^ and q[" 4BC '' ) are the calcula- 
tions with reflecting and absorbing boundaries respect- 
ively. This is to see how any error from the absorbing 
boundary effects the interior points. Figure 2 shows the 
result for each angular momentum with two different grid 
spacings. Again the reflecting boundaries are chosen to 
be at r = 100. 

We see that in all cases the error has remained small 
throughout the interior, for the dx — dt — 0.1 case 
bounded by 10 -3 and for dx = dt = 0.01 bounded by 
10~ 5 . This is within the 0(Ar 2 , At 2 ) expected from the 
discretisation. 



C. Temporal Comparison of Probability 

We now test the how the error evolves through time. 
This is done by calculating the probability of finding the 
particle inside the interval over time, mathematically the 
following is calculated: 



P{t) 



10 



\Qi(r,t)\ 2 dr 



(62) 



with reflecting and absorbing boundaries and the abso- 
lute value of the difference taken. 

For this test we increase the time interval to [0, 50] and 
move the reflecting boundary to r = 200. In each case 
more than 90% of the wavefunction has left the inter- 
val, specifically the probabilities inside the interval are 
8.57E - 002, 6.36E - 003 and 2.03E - 004 for I = 0, 1, 2 
respectively at the end of the calculation. 

Figure 3 shows the results for each angular momenta 
and different grid spacings. 

We see that in time also the error remains bounded. 
From the plots it appears the bound on the error is pro- 
portional to the grid spacings. 

These results are satisfactory and so now with confid- 
ence in the previous work we go on to the Hartree-Fock 
calculations. 
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Figure 2. The figures shows a comparison of the radial com- 
ponent of the wavefunctions at the final time 15, for angular 
momenta I = 0,1,2, calculated with each technique. The 
value in equation (61) is plotted against the radius. 



VIII. RESULTS AND TESTING: 
HARTREE-FOCK RESONANCES IN THE 
LINEAR REGIME 
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Results from the implementation of the discretised 
Hartree-Fock system, as described in sections IV and VI, 
are now shown. We first present the variation of the root 
mean square radius over time for ^He, 1 ^0 and For 
each nuclei the following is shown: 

(a) A calculation performed with reflecting boundaries 
at 1500 ha. This is the result expected from a 
continuum calculation because the boundary is far 
enough away so as to avoid reflection. This is plot- 
ted from to 500 fm c _1 to show the main features 
occurring at the beginning of the resonance. 



Figure 3. (Color online) These plots show how the error in the 
probability from the absorbing boundaries changes through 
time. Equation (62) is calculated with reflecting and ab- 
sorbing boundaries and the absolute value of there difference 
taken, though time and plotted. 



(b) The result of using reflecting boundaries at 30 fm. 
This is to show the effect the absorbing boundar- 
ies are having. Again this is plotted from to 
500 fm c" 1 . 
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(c) The difference between the expected result in (a) and 
a calculation with absorbing boundaries at 30 fm. 
This is plotted for the entire to 3000 fm c _1 time 
range. This difference is an error due to the discretiz- 
ation of the absorbing boundaries and so we consider 
a upper bound for this value of 0(Ar 2 ) acceptable. 

For each nucleus has a group of three figures are shown 
which are labelled according to the above. We also show 
the time each calculation takes to evaluate the efficiency 
of the absorbing bounds. 

Grid spacings of Ar = 0.1 fm and At = 0.1 fm c _1 are 
used and all calculation are evolved from to 3000 fm 



1 . Hehum-4 

From figure (4a) we can see that the resonance for ^He 
has a simple damped oscillatory motion, the radius of 
the nuclei repeatedly increasing and decreasing clearly 
demonstrating the breathing mode. Figure (4c) shows us 
that the absorbing boundary provide us with a reasonable 
discrepancy from the expected result being bounded by 
10~ 7 , well below the O(0. 1 2 ) discretization error. Finally 
by comparing (4a) and (4b) the effect of the reflected flux 
can clearly be seen, which is the source of discretisation 
artefacts in the strength functions [44]. 



2. Oxygen-16 

The top panel of figure (5a) shows a more complicated 
motion of the nucleus this time, which does not look like 
a single damped mode. This is due to the multiple single- 
particle states present, known as Landau fragmentation. 
The absolute error as shown in figure (5c) is bounded 
by a larger number than helium, but again within the 
acceptable range. 



3. Calcium- 40 

The results for calcium again show a damped oscilla- 
tion, as expected, though a long-lived resonant compon- 
ent is excited too, which the reflecting boundaries ob- 
viously cannot reproduce for long times. The errors are 
somewhat larger than the helium or oxygen cases but still 
acceptable. 



A. Timing 

As an guide, we present a table of timing results for 
the Oxygen calculations in Table I. 

The results show that the absorbing boundaries are 
considerably more expensive than reflecting boundaries, 
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Figure 4. The time evolution of the monopole moment in 
Helium-4, showing (a) the continuum result, (b) for compar- 
ison, the result of a reflecting boundary wall and (c) the ab- 
solute value of the difference between the monopole moments 
when calculated using a absorbing boundary and using a far 
reflecting wall, over time. 



but less so than using a large box with simple bound- 
ary conditions. It is interesting also to examine the time 
taken to each iteration. Figure (7) shows a plot of the 
time to compute each iteration, as a running average over 
20 iterations to somewhat smooth out the effect of com- 
puter load. This shows the steady increase in expense to 
calculate a iteration as the calculation progresses and is 
due to the non-locality in time of the absorbing boundary 
condition. 
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Figure 5. The time evolution of the monopole moment in 
Oxygen-16, showing (a) the continuum result, (b) for com- 
parison, the result of a reflecting boundary wall and (c) the 
absolute value of the difference between the monopole mo- 
ments when calculated using a absorbing boundary and using 
a far reflecting wall, over time. 



Figure 6. The time evolution of the monopole moment in 
Calcium-40, showing (a) the continuum result, (b) for com- 
parison, the result of a reflecting boundary wall and (c) the 
absolute value of the difference between the monopole mo- 
ments when calculated using a absorbing boundary and using 
a far reflecting wall, over time. 



B. Strength Functions 

The strength functions for these calculations are now 
presented. As these are the calculations required in order 
to make comparisons to experiment their accurate calcu- 
lation is critical. We require that the error in the above 
results do not give noticable artefacts in the strength 
functions, at least to the level of experimental resolu- 
tion. Figure (8) shows the calculated strength function 
from the expected result with that calculated using ab- 
sorbing boundaries. 

We see that both calculations match up well for all the 



nuclei tested. The figures show the increasing complexity 
of the nuclear structure, as more features appear in the 
strength functions. 



IX. RESULTS AND TESTING: NON-LINEAR 
REGIME 

As well as testing in the small amplitude linear re- 
sponse regime, of relevance to giant resonances, it is 
also instructive to examine the larger-amplitude regime, 
which can be studied in THDF-based techniques [45-47] , 
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Boundary Type 


R(fm) 


Calculation Time (s) 


Reflecting 


1500 


2378 


Reflecting 


30 


58 


Absorbing 


30 


144 



Table I. Calculation times for the large box continuum cal- 
culation with reflecting bounds, a small-box calculation with 
spurious reflections and a small-box calculation with absorb- 
ing boundaries. 
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Figure 7. A plot showing the expense of each iteration in a 
calculation of oxygen-16. It clearly show the non-locally of 
the absorbing boundary increasing the calculation time for 
iteration the further the calculation progresses. 



unlike the small-amplitude-limited RPA. This regime is 
relevant to the decay of highly excited fragments follow- 
ing e.g. deep inelastic collisions, and significant particle 
emission may be expected. Similar situation arise in 
atomic physics where direct electromagnetic excitation of 
highly ionizing collective modes is feasible [48]. We use 
a test case of monopole exciations of 16 O, with increas- 
ingly strong boosts (9) such that eventually all particles 
are lost from the nucleus through large-amplitude excit- 
ation. We note that the computational effort for large 
amplitude excitations is not different to that for small- 
amplitude excitations, as the iteration procedure is not 
changed for larger amplitudes. 

Despite the success of the small-amplitude calcula- 
tions, there is no a priori reason to expect larger amp- 
litude calculations to perform so well, since our absorbing 
boundaries are predicated on the fact that the only po- 
tential active at the boundary is the centrifugal barrier, 
whereas the nuclear mean-field exists wherever the nuc- 
leon wavefunction is finite. As more particles are emitted, 
so too the nuclear wavefunction and its associated mean- 
field are present in the exterior region. Figure 9 shows 
the comparison of the total number of particles emitted 
(by 1500 fm/c) from a 16 nucleus between an absorbing 
boundary calculation, and a reflecting boundary calcula- 
tion in which the size of the box is so high that the re- 
flecting boundaries are not reached. The range of boost 
is sufficiently large to cover the small amplitude limit as 
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Figure 8. Plots showing the effect of using the absorbing 
boundary condition on the strength functions of various nuc- 
lei. Going from top to bottom there is the helium, oxygen 
and calium strength functions. 



well as the regime in which the nucleus is entirely ionized. 
The two calculations are seen to be close over the entire 
range, with small differences near the bend as complete 
ionization occurs. The time-dependence of the particle 
emission is shown in Figure 10, in which the case around 
the bend is shown to still be changing at the end time of 
the calculation. 

Figure 11 shows the time-dependent error (absorbing 
bounds compared with large-space reflecting bounds) in 
the total number of particles emitted for a range of kick 
size. This highlights the small differences in Figure 9 
where the errors around k = 0.2 fm~ 2 are seen to be 
largest. In the worst case, this error is noticable, but 
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Figure 9. A comparison of the number of particles emitted 
from the region between and 30fm with absorbing bound- 
aries at 30fm compared with reflecting boundaries at 600fm 
which are not reached in the time of the calculation. 
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Figure 10. The time-dependence of particle emission as a 
function of boost strength for large-amplitude excitations in 
16 0. The legend indicates the strength k (fin -2 ) of the boost 
in equation (9). 



still rather small. 



X. PERSPECTIVES AND CONCLUSION 



A. Perspectives for more realistic calculations 



Our calculations respresent a step on the way to more 
realistic calculations of giant resonances within a con- 
tinuum time-dependent Hartree-Fock framework. We 
discuss in this section some perspectives for the possibil- 
ity of performing more realistic calculations. Our calcu- 
lations deliberately considered a simple case, yet within 
TDHF-based methods, calculations without our form of 
absorbing bound exist with more relaxed symmetries 
[38, 49-51] or with pairing in the BCS or TDHFB frame- 
work [49, 52-54]. Our method is extendable in a straight- 
forward way to calculations involving pairing. The in- 
creased expense scales in the same way as discrete cal- 
culations with pairing scale with respect to calculations 
without pairing. The addition of extra single-particle 
states to account for the scattering of Cooper pairs will 
involve extra boundary conditions, but only with a linear 
scaling with respect to the number of particle states. On 
the other hand, increased dimensions will be more costly. 
In our case of spherical symmetry in which there is a 
single boundary point for 300 interior points, we have a 
similar time spent on the boundary as the entire internal 
region. In a three-dimensional calculation, in which the 
boundary is the surface of volume, the ratio of boundary 
points to internal points is much higher. Our technique 
is thus not currently suitable for a three-dimensional cal- 
culation. However, reasonable scaling could nevertheless 
be achieved with an expansion of the density in spherical 
harmonics. For the purposes of calculating giant reson- 
ances of general multipolarity and of deformed nuclei, 
this would suffice, as only one point per moment of the 
density would be needed to act as a boundary point, and 
a typical expansion of a handful of terms would describe 
a small-amplitude deformation. A full three-dimensional 
code would remain required for heavy-ion collisions. 

Our immediate aim is to find a suitable way to in- 
clude the Coulomb potential, which has been ignored 
here, within the treatment of the absorbing boundaries. 
The practical realisation of this is more difficult than the 
present case because the required inverse Laplace trans- 
form is not of a simple form. The current approach being 
developed is to use the method in [55-57] to approximate 
the more complex inverse Laplace transform. 

It should also be possible to reduce the time taken to 
perform the boundary calculation. In the oxygen tests it 
was shown that most of the expense comes from the end 
of the calculation where the non-locality in time plays a 
part. One solution to this would be to use the method de- 
scribed in [58] which uses a sum of exponentials approx- 
imation that can be evaluated recursively. The effect is 
to reduce the sum in (57) that requires O(N) operations 
to one that requires just 0(lniV). 
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Figure 11. The total error in number of particles emitted by the nucleus as a function of time for increasingly stronger boosts 
(indicated by the strength k in each panel). The error is calculated with respect to a calculation withouth aborsbing bounds 
but in a space so large that the boundaries are not probed. The boost parameter k is as defined in (9). 



B. Conclusion 

We have presented a implementation of a spheric- 
ally symmetric Hartree-Fock system discretised using a 
Crank-Nicholson scheme. We also presented the deriva- 
tion and implementation of an absorbing boundary con- 
dition approach to handle the outgoing wave condition. 
It was shown using a Laplace transform method that it 
is possible to construct a boundary condition at a finite 
distance away from the origin. This came at the cost of it 
being non-local in time, meaning the value of the wave- 
function at the boundary has to be stored throughout 
the calculation, causing an increase in the time taken to 
calculate each iteration as it progressed. 



The results of the testing show that absorbing bound- 
ary conditions do provide a suitable way of treating the 
boundary in spatially unbounded time-dependent prob- 
lems. We see that although there are errors introduced 
from the discretization of the absorbing boundaries, they 
are small and stay small throughout the various manip- 
ulations required to calculate the strength functions. As 
well as being accurate they also show a good improve- 
ment in the speed of the calculation compared to using 
a large box. 

We applied the method to large amplitude motion, and 
found acceptable results. We discussed perspectives for 
future, and more realistic, calculations. 
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